suppressWarnings(library(MendelianRandomization))
suppressWarnings(library(digest))
suppressWarnings(library(tidyr))
suppressWarnings(library(MRPRESSO))
suppressWarnings(library(rowr))
suppressWarnings(library(ggplot2))
suppressWarnings(library(dplyr))
setwd('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants')
diabetes_all<-read.csv('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Diabetes/Newrf_diabetes_instruments.csv', header = TRUE, stringsAsFactors = FALSE)
#diabetes_all<-diabetes_all[diabetes_all$F.statistic>10,] Use command to remove weak instruments
breast<-readRDS('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Breast Cancer/breast.rds')
prostate<-readRDS('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Prostate Cancer/prostate.rds')
breast= separate(data=breast, col=phase3_1kg_id, into=c('SNP','rest'), sep= "\\:") #Create SNP column in breast data
#Merge diabetes and breast data
index_b<-match(breast$SNP,diabetes_all$SNP..rs.id.)
diabetes_b<-cbind(diabetes_all[index_b,],breast)
#ALign effect alleles in 2 datasets
data_in_ea = diabetes_b$Effect.Allele
data_in_nea = diabetes_b$Other.Allele
ea_new = diabetes_b$a1
nea_new = diabetes_b$a0
stay=as.numeric(data_in_ea == ea_new) #stay = 1 if the same
swap=as.numeric(data_in_nea == ea_new)
#checks
table(stay)
## stay
## 0 1
## 57 54
table(swap)
## swap
## 0 1
## 54 57
stay*swap == 0
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [99] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
swap_vec = rep(1, length(swap))
swap_vec = (-2 * swap) + swap_vec
#align the beta vector in the new data
diabetes_b$bcac_onco_icogs_gwas_beta =swap_vec*as.numeric(diabetes_b$bcac_onco_icogs_gwas_beta)
diabetes_b$Effect = swap_vec * as.numeric(diabetes_b$Effect)
####################################################################################
#Merge diabetes and prostate data
index_p<-match(prostate$SNP,diabetes_all$SNP..rs.id.)
diabetes_p<-cbind(diabetes_all[index_p,],prostate)
colnames(diabetes_p)[23]<-'Effect_p'
colnames(diabetes_p)[24]<-'StdErr_p'
#ALign effect alleles in 2 datasets
data_in_ea = diabetes_p$Effect.Allele
data_in_nea = diabetes_p$Other.Allele
ea_new = diabetes_p$Allele1
nea_new = diabetes_p$Allele2
stay=as.numeric(data_in_ea == ea_new) #stay = 1 if the same
swap=as.numeric(data_in_nea == ea_new)
#checks
table(stay)
## stay
## 0
## 88
table(swap)
## swap
## 0
## 88
stay*swap == 0
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [85] TRUE TRUE TRUE TRUE
swap_vec = rep(1, length(swap))
swap_vec = (-2 * swap) + swap_vec
#align the beta vector in the new data
diabetes_p$Effect_p =swap_vec*as.numeric(diabetes_p$Effect_p)
diabetes_p$Effect = swap_vec * as.numeric(diabetes_p$Effect)
# remove NAs, we believe that they are missing at random
diabetes_b<-na.omit(diabetes_b)
diabetes_p<-na.omit(diabetes_p)
#Get rid of duplicated columns
diabetes_b <- diabetes_b[, !duplicated(colnames(diabetes_b))]
diabetes_p <- diabetes_p[, !duplicated(colnames(diabetes_p))]
#Load data of rs ids of grouped SNPs from PhenoScanner query
BMI=read.csv('C:/Users/marta/Desktop/IC/Thesis/BMI.csv', stringsAsFactors = FALSE, header=FALSE)
Glucose=read.csv('C:/Users/marta/Desktop/IC/Thesis/Glucose.csv',stringsAsFactors = FALSE,header=FALSE)
Insulin=read.csv('C:/Users/marta/Desktop/IC/Thesis/Insulin.csv',stringsAsFactors = FALSE, header=FALSE)
Lipids=read.csv('C:/Users/marta/Desktop/IC/Thesis/Lipids.csv',stringsAsFactors = FALSE, header=FALSE)
MR breast cancer
MRInputObject_breast <- mr_input(bx = diabetes_b$Effect,
bxse = diabetes_b$StdErr,
by = diabetes_b$bcac_onco_icogs_gwas_beta,
byse = diabetes_b$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=diabetes_b$SNP)
MRAllObject_all_breast <- mr_allmethods(MRInputObject_breast, method = "all")
MRAllObject_all_breast
## Method Estimate Std Error 95% CI P-value
## Simple median -0.004 0.019 -0.041 0.033 0.835
## Weighted median 0.011 0.018 -0.024 0.045 0.539
## Penalized weighted median -0.004 0.018 -0.039 0.031 0.825
##
## IVW 0.000 0.017 -0.033 0.033 0.998
## Penalized IVW -0.007 0.013 -0.033 0.019 0.581
## Robust IVW 0.005 0.023 -0.040 0.050 0.819
## Penalized robust IVW -0.008 0.016 -0.040 0.024 0.612
##
## MR-Egger 0.029 0.035 -0.040 0.098 0.407
## (intercept) -0.003 0.003 -0.009 0.003 0.344
## Penalized MR-Egger 0.003 0.031 -0.058 0.064 0.926
## (intercept) -0.001 0.002 -0.006 0.004 0.743
## Robust MR-Egger 0.046 0.066 -0.083 0.176 0.484
## (intercept) -0.004 0.005 -0.013 0.006 0.447
## Penalized robust MR-Egger -0.005 0.041 -0.085 0.076 0.908
## (intercept) 0.000 0.003 -0.006 0.006 0.939
mr_plot(MRInputObject_breast,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast)
#MR-PRESSO breast
input = data.frame(betaY=diabetes_b$bcac_onco_icogs_gwas_beta, seY=diabetes_b$bcac_onco_icogs_gwas_se,beta_diabetes=diabetes_b$Effect, se_diabetes=diabetes_b$StdErr,row.names= NULL)
mrpresso_out = mr_presso(BetaOutcome = "betaY", BetaExposure = c("beta_diabetes"), SdOutcome = "seY", SdExposure="se_diabetes", data=input, OUTLIERtest = TRUE, DISTORTIONtest = TRUE, NbDistribution = 1000)
## Warning in mr_presso(BetaOutcome = "betaY", BetaExposure =
## c("beta_diabetes"), : Outlier test unstable. The significance threshold
## of 0.05 for the outlier test is not achievable with only 1000 to
## compute the null distribution. The current precision is <0.111. Increase
## NbDistribution.
mrpresso_out$`MR-PRESSO results`$`Global Test`$`RSSobs`
## [1] 459.3003
mrpresso_out$`MR-PRESSO results`$`Global Test`$Pvalue
## [1] "<0.001"
mrpresso_out$`Main MR results`
## Exposure MR Analysis Causal Estimate Sd T-stat
## 1 beta_diabetes Raw 4.022679e-05 0.01704120 0.002360562
## 2 beta_diabetes Outlier-corrected -7.200820e-03 0.01400442 -0.514181786
## P-value
## 1 0.9981208
## 2 0.6082263
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Distortion Coefficient`
## beta_diabetes
## 100.5586
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$Pvalue
## [1] 0.228
out_indices = mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Outliers Indices`
outliers = diabetes_b$SNP[out_indices]
# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%outliers,yes='red',no='blue'))) +
geom_point(size=1.5) +
labs(title = "", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "outliers"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 10), axis.title.x = element_text(size = 12),
axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 12),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%BMI$V2,yes='red',no='blue'))) +
geom_point(size=1.5) +
labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "BMI SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
index_bmi<-match(diabetes_b$SNP,BMI$V2)
data_b_bmi<-na.omit(diabetes_b[index_bmi,])
MRInputObject_breast_bmi <- mr_input(bx = data_b_bmi$Effect,
bxse = data_b_bmi$StdErr,
by = data_b_bmi$bcac_onco_icogs_gwas_beta,
byse = data_b_bmi$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_bmi$SNP)
MRAllObject_all_breast_bmi<- mr_allmethods(MRInputObject_breast_bmi, method = "all")
MRAllObject_all_breast_bmi
## Method Estimate Std Error 95% CI P-value
## Simple median -0.034 0.029 -0.091 0.023 0.242
## Weighted median -0.018 0.027 -0.070 0.035 0.508
## Penalized weighted median -0.028 0.027 -0.082 0.025 0.299
##
## IVW -0.010 0.025 -0.059 0.039 0.696
## Penalized IVW -0.031 0.021 -0.072 0.010 0.138
## Robust IVW -0.010 0.031 -0.071 0.050 0.739
## Penalized robust IVW -0.031 0.025 -0.080 0.019 0.224
##
## MR-Egger 0.052 0.055 -0.056 0.161 0.344
## (intercept) -0.006 0.004 -0.014 0.003 0.209
## Penalized MR-Egger 0.064 0.048 -0.031 0.158 0.188
## (intercept) -0.006 0.004 -0.014 0.002 0.128
## Robust MR-Egger 0.058 0.065 -0.070 0.186 0.375
## (intercept) -0.006 0.005 -0.016 0.004 0.215
## Penalized robust MR-Egger 0.066 0.055 -0.042 0.174 0.233
## (intercept) -0.006 0.004 -0.014 0.002 0.134
mr_plot(MRInputObject_breast_bmi,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_bmi)
# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Insulin$V2,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Insulin SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
index<-match(diabetes_b$SNP,Insulin$V2)
data_b_insulin<-na.omit(diabetes_b[index,])
MRInputObject_breast_insulin <- mr_input(bx = data_b_insulin$Effect,
bxse = data_b_insulin$StdErr,
by = data_b_insulin$bcac_onco_icogs_gwas_beta,
byse = data_b_insulin$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_insulin$SNP)
MRAllObject_all_breast_insulin <- mr_allmethods(MRInputObject_breast_insulin, method = "all")
MRAllObject_all_breast_insulin
## Method Estimate Std Error 95% CI P-value
## Simple median -0.072 0.069 -0.207 0.063 0.297
## Weighted median 0.082 0.047 -0.010 0.175 0.081
## Penalized weighted median 0.132 0.047 0.040 0.224 0.005
##
## IVW 0.026 0.062 -0.095 0.146 0.674
## Penalized IVW -0.043 0.047 -0.135 0.048 0.350
## Robust IVW 0.027 0.061 -0.093 0.146 0.663
## Penalized robust IVW -0.041 0.064 -0.167 0.085 0.520
##
## MR-Egger 0.203 0.111 -0.015 0.422 0.068
## (intercept) -0.014 0.008 -0.029 0.001 0.064
## Penalized MR-Egger 0.284 0.076 0.134 0.434 0.000
## (intercept) -0.020 0.005 -0.030 -0.010 0.000
## Robust MR-Egger 0.221 0.104 0.016 0.425 0.034
## (intercept) -0.016 0.007 -0.030 -0.002 0.027
## Penalized robust MR-Egger 0.286 0.080 0.129 0.443 0.000
## (intercept) -0.020 0.006 -0.032 -0.008 0.001
mr_plot(MRInputObject_breast_insulin,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_insulin)
# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Glucose$V2,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Glucose SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
index<-match(diabetes_b$SNP,Glucose$V2)
data_b_glucose<-na.omit(diabetes_b[index,])
MRInputObject_breast_glucose <- mr_input(bx = data_b_glucose$Effect,
bxse = data_b_glucose$StdErr,
by = data_b_glucose$bcac_onco_icogs_gwas_beta,
byse = data_b_glucose$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_glucose$SNP)
MRAllObject_all_breast_glucose <- mr_allmethods(MRInputObject_breast_glucose, method = "all")
MRAllObject_all_breast_glucose
## Method Estimate Std Error 95% CI P-value
## Simple median -0.072 0.076 -0.221 0.078 0.346
## Weighted median -0.001 0.058 -0.114 0.113 0.989
## Penalized weighted median -0.070 0.061 -0.189 0.049 0.248
##
## IVW -0.004 0.071 -0.143 0.136 0.960
## Penalized IVW -0.078 0.045 -0.167 0.011 0.087
## Robust IVW -0.010 0.080 -0.168 0.147 0.897
## Penalized robust IVW -0.084 0.068 -0.218 0.050 0.220
##
## MR-Egger 0.067 0.135 -0.198 0.331 0.621
## (intercept) -0.005 0.008 -0.022 0.011 0.536
## Penalized MR-Egger 0.025 0.097 -0.165 0.216 0.794
## (intercept) -0.003 0.006 -0.016 0.009 0.610
## Robust MR-Egger 0.063 0.135 -0.203 0.328 0.642
## (intercept) -0.005 0.009 -0.023 0.012 0.539
## Penalized robust MR-Egger 0.009 0.162 -0.309 0.327 0.955
## (intercept) -0.003 0.008 -0.019 0.012 0.690
mr_plot(MRInputObject_breast_glucose,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_glucose)
index<-match(diabetes_b$SNP,Lipids$V2)
data_b_lipid<-na.omit(diabetes_b[index,])
# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Lipids$V2,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Lipid SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
MRInputObject_breast_lipid <- mr_input(bx = data_b_lipid$Effect,
bxse = data_b_lipid$StdErr,
by = data_b_lipid$bcac_onco_icogs_gwas_beta,
byse = data_b_lipid$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_lipid$SNP)
MRAllObject_all_breast_lipid <- mr_allmethods(MRInputObject_breast_lipid, method = "all")
MRAllObject_all_breast_lipid
## Method Estimate Std Error 95% CI P-value
## Simple median -0.077 0.089 -0.252 0.099 0.390
## Weighted median -0.075 0.066 -0.205 0.055 0.257
## Penalized weighted median -0.130 0.065 -0.257 -0.003 0.045
##
## IVW -0.011 0.080 -0.168 0.146 0.894
## Penalized IVW -0.105 0.049 -0.202 -0.009 0.033
## Robust IVW -0.025 0.109 -0.239 0.189 0.820
## Penalized robust IVW -0.130 0.064 -0.256 -0.005 0.042
##
## MR-Egger 0.059 0.144 -0.223 0.342 0.680
## (intercept) -0.005 0.009 -0.023 0.012 0.554
## Penalized MR-Egger 0.012 0.104 -0.192 0.217 0.905
## (intercept) -0.004 0.007 -0.017 0.010 0.595
## Robust MR-Egger 0.052 0.150 -0.242 0.345 0.731
## (intercept) -0.006 0.009 -0.023 0.012 0.541
## Penalized robust MR-Egger -0.053 0.220 -0.484 0.379 0.811
## (intercept) -0.002 0.007 -0.017 0.012 0.742
mr_plot(MRInputObject_breast_lipid,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_lipid)
BMI=read.csv('C:/Users/marta/Desktop/IC/Thesis/Genetic variants/BMI.csv', stringsAsFactors = FALSE, header=FALSE)
index<-match(diabetes_b$SNP,BMI$V1)
data_b_bmi<-cbind(diabetes_b,BMI[index,])
data_b_bmi$V2[is.na(data_b_bmi$V2)] <- 'NO BMI'
data_nobmi<-data_b_bmi[data_b_bmi$V2=='NO BMI',]
MRInputObject_breast_nobmi <- mr_input(bx = data_nobmi$Effect,
bxse = data_nobmi$StdErr,
by = data_nobmi$bcac_onco_icogs_gwas_beta,
byse = data_nobmi$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_nobmi$SNP)
MRAllObject_all_breast_nobmi <- mr_allmethods(MRInputObject_breast_nobmi, method = "all")
MRAllObject_all_breast_nobmi
## Method Estimate Std Error 95% CI P-value
## Simple median 0.036 0.030 -0.022 0.095 0.221
## Weighted median -0.005 0.027 -0.059 0.049 0.856
## Penalized weighted median -0.004 0.028 -0.058 0.050 0.871
##
## IVW -0.012 0.022 -0.056 0.031 0.583
## Penalized IVW 0.002 0.021 -0.038 0.043 0.906
## Robust IVW -0.010 0.022 -0.052 0.033 0.649
## Penalized robust IVW 0.003 0.020 -0.036 0.041 0.898
##
## MR-Egger -0.040 0.048 -0.134 0.055 0.408
## (intercept) 0.002 0.003 -0.004 0.009 0.517
## Penalized MR-Egger -0.043 0.046 -0.133 0.048 0.358
## (intercept) 0.003 0.003 -0.004 0.009 0.374
## Robust MR-Egger -0.047 0.041 -0.127 0.034 0.254
## (intercept) 0.003 0.004 -0.004 0.011 0.400
## Penalized robust MR-Egger -0.049 0.039 -0.126 0.029 0.217
## (intercept) 0.004 0.003 -0.003 0.010 0.296
mr_plot(MRInputObject_breast_nobmi,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_nobmi)
MRInputObject_prostate <- mr_input(bx = diabetes_p$Effect,
bxse = diabetes_p$StdErr,
by = diabetes_p$Effect_p,
byse = diabetes_p$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=diabetes_p$SNP)
MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate, method = "all")
MRAllObject_all_prostate
## Method Estimate Std Error 95% CI P-value
## Simple median 0.011 0.024 -0.035 0.058 0.632
## Weighted median 0.001 0.022 -0.042 0.045 0.947
## Penalized weighted median -0.005 0.022 -0.049 0.039 0.819
##
## IVW -0.041 0.042 -0.124 0.042 0.333
## Penalized IVW -0.010 0.016 -0.042 0.022 0.543
## Robust IVW -0.007 0.018 -0.042 0.028 0.690
## Penalized robust IVW -0.012 0.016 -0.043 0.020 0.469
##
## MR-Egger -0.113 0.100 -0.309 0.084 0.260
## (intercept) 0.006 0.007 -0.009 0.021 0.428
## Penalized MR-Egger -0.048 0.038 -0.122 0.027 0.213
## (intercept) 0.003 0.003 -0.003 0.009 0.297
## Robust MR-Egger -0.020 0.039 -0.097 0.056 0.603
## (intercept) 0.001 0.003 -0.004 0.006 0.690
## Penalized robust MR-Egger -0.045 0.032 -0.108 0.017 0.155
## (intercept) 0.003 0.002 -0.002 0.007 0.274
mr_plot(MRInputObject_prostate,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)
#MR-PRESSO prostate
input = data.frame(betaY=diabetes_p$Effect_p,seY=diabetes_p$StdErr_p,beta_diabetes=diabetes_p$Effect,se_diabetes=diabetes_p$StdErr, row.names= NULL)
mrpresso_out = mr_presso(BetaOutcome = "betaY", BetaExposure = "beta_diabetes", SdOutcome = "seY", SdExposure="se_diabetes", data=input, OUTLIERtest = TRUE, DISTORTIONtest = TRUE, NbDistribution = 1000)
## Warning in mr_presso(BetaOutcome = "betaY", BetaExposure =
## "beta_diabetes", : Outlier test unstable. The significance threshold of
## 0.05 for the outlier test is not achievable with only 1000 to compute
## the null distribution. The current precision is <0.088. Increase
## NbDistribution.
mrpresso_out$`Main MR results`
## Exposure MR Analysis Causal Estimate Sd T-stat
## 1 beta_diabetes Raw -0.04085910 0.04223478 -0.9674276
## 2 beta_diabetes Outlier-corrected -0.01429709 0.01722820 -0.8298654
## P-value
## 1 0.3360120
## 2 0.4089954
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Distortion Coefficient`
## beta_diabetes
## -185.7861
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$Pvalue
## [1] 0.044
out_indices = mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Outliers Indices`
outliers=diabetes_p$SNP[out_indices]
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%outliers,yes='red',no='blue'))) +
geom_point(size=1.5) +
labs(title='',x = "Genetic association with T2D", y = "Genetic association with prostate cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "outliers"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 10), axis.title.x = element_text(size = 12),
axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 12),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
BMI=read.csv('C:/Users/marta/Desktop/IC/Thesis/BMI.csv', stringsAsFactors = FALSE, header=FALSE)
index<-match(diabetes_p$SNP,BMI$V2)
data_p_bmi<-na.omit(diabetes_p[index,])
Creating scatter plot
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%BMI$V2,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "BMI SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
MRInputObject_prostate_bmi <- mr_input(bx = data_p_bmi$Effect,
bxse = data_p_bmi$StdErr,
by = data_p_bmi$Effect_p,
byse = data_p_bmi$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_bmi$SNP)
MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_bmi, method = "all")
MRAllObject_all_prostate
## Method Estimate Std Error 95% CI P-value
## Simple median -0.019 0.034 -0.085 0.047 0.569
## Weighted median -0.036 0.032 -0.099 0.026 0.257
## Penalized weighted median -0.048 0.032 -0.111 0.015 0.135
##
## IVW -0.016 0.027 -0.070 0.038 0.561
## Penalized IVW -0.032 0.024 -0.078 0.015 0.179
## Robust IVW -0.026 0.026 -0.077 0.025 0.322
## Penalized robust IVW -0.034 0.023 -0.079 0.011 0.139
##
## MR-Egger -0.006 0.073 -0.150 0.138 0.934
## (intercept) -0.001 0.006 -0.013 0.011 0.884
## Penalized MR-Egger -0.050 0.064 -0.175 0.075 0.433
## (intercept) 0.002 0.005 -0.009 0.012 0.759
## Robust MR-Egger -0.027 0.053 -0.132 0.078 0.611
## (intercept) 0.000 0.004 -0.008 0.008 0.961
## Penalized robust MR-Egger -0.049 0.042 -0.132 0.034 0.244
## (intercept) 0.001 0.004 -0.006 0.009 0.701
mr_plot(MRInputObject_prostate_bmi,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)
index<-match(diabetes_p$SNP,Insulin$V2)
data_p_insulin<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Insulin$V2,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Insulin SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
MRInputObject_prostate_insulin <- mr_input(bx = data_p_insulin$Effect,
bxse = data_p_insulin$StdErr,
by = data_p_insulin$Effect_p,
byse = data_p_insulin$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_insulin$SNP)
MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_insulin, method = "all")
## Warning in lmrob.fit(x, y, control, init = init): M-step did NOT converge.
## Returning unconverged SM-estimate
MRAllObject_all_prostate
## Method Estimate Std Error 95% CI P-value
## Simple median 0.084 0.096 -0.103 0.272 0.378
## Weighted median 0.088 0.078 -0.066 0.242 0.262
## Penalized weighted median 0.088 0.083 -0.074 0.250 0.288
##
## IVW 0.078 0.081 -0.080 0.236 0.333
## Penalized IVW 0.098 0.060 -0.019 0.215 0.100
## Robust IVW 0.091 0.232 -0.364 0.547 0.694
## Penalized robust IVW 0.090 0.064 -0.035 0.216 0.159
##
## MR-Egger 0.069 0.162 -0.249 0.388 0.669
## (intercept) 0.001 0.011 -0.021 0.022 0.949
## Penalized MR-Egger 0.068 0.122 -0.170 0.306 0.577
## (intercept) 0.001 0.007 -0.013 0.014 0.914
## Robust MR-Egger 0.048 NA NA NA NA
## (intercept) 0.002 NA NA NA NA
## Penalized robust MR-Egger 0.066 0.126 -0.180 0.312 0.599
## (intercept) 0.001 0.005 -0.009 0.011 0.870
mr_plot(MRInputObject_prostate_insulin,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)
index<-match(diabetes_p$SNP,Glucose$V2)
data_p_glucose<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Glucose$V2,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Glucose SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
MRInputObject_prostate_glucose <- mr_input(bx = data_p_glucose$Effect,
bxse = data_p_glucose$StdErr,
by = data_p_glucose$Effect_p,
byse = data_p_glucose$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_glucose$SNP)
MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_glucose, method = "all")
MRAllObject_all_prostate
## Method Estimate Std Error 95% CI P-value
## Simple median 0.060 0.107 -0.148 0.269 0.571
## Weighted median 0.065 0.080 -0.093 0.222 0.420
## Penalized weighted median 0.027 0.077 -0.125 0.179 0.727
##
## IVW 0.135 0.064 0.009 0.261 0.036
## Penalized IVW 0.067 0.060 -0.051 0.185 0.266
## Robust IVW 0.036 0.063 -0.087 0.160 0.562
## Penalized robust IVW 0.055 0.062 -0.066 0.176 0.376
##
## MR-Egger 0.238 0.110 0.023 0.453 0.030
## (intercept) -0.007 0.006 -0.020 0.005 0.252
## Penalized MR-Egger 0.238 0.110 0.023 0.453 0.030
## (intercept) -0.007 0.006 -0.020 0.005 0.252
## Robust MR-Egger 0.167 0.411 -0.640 0.973 0.685
## (intercept) -0.005 0.014 -0.033 0.023 0.732
## Penalized robust MR-Egger 0.167 0.411 -0.640 0.973 0.685
## (intercept) -0.005 0.014 -0.033 0.023 0.732
mr_plot(MRInputObject_prostate_glucose,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)
index<-match(diabetes_p$SNP,Lipids$V2)
data_p_lipid<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Lipids$V2,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Lipid SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
MRInputObject_prostate_lipid <- mr_input(bx = data_p_lipid$Effect,
bxse = data_p_lipid$StdErr,
by = data_p_lipid$Effect_p,
byse = data_p_lipid$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_lipid$SNP)
MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_lipid, method = "all")
MRAllObject_all_prostate
## Method Estimate Std Error 95% CI P-value
## Simple median 0.075 0.103 -0.127 0.277 0.465
## Weighted median 0.083 0.093 -0.099 0.266 0.372
## Penalized weighted median 0.053 0.090 -0.124 0.229 0.558
##
## IVW 0.154 0.074 0.008 0.300 0.038
## Penalized IVW 0.090 0.067 -0.042 0.222 0.182
## Robust IVW 0.114 0.173 -0.225 0.453 0.508
## Penalized robust IVW 0.081 0.068 -0.053 0.215 0.238
##
## MR-Egger 0.417 0.120 0.182 0.653 0.001
## (intercept) -0.019 0.008 -0.034 -0.004 0.012
## Penalized MR-Egger 0.417 0.120 0.182 0.653 0.001
## (intercept) -0.019 0.008 -0.034 -0.004 0.012
## Robust MR-Egger 0.420 0.122 0.181 0.660 0.001
## (intercept) -0.019 0.007 -0.033 -0.006 0.006
## Penalized robust MR-Egger 0.420 0.122 0.181 0.660 0.001
## (intercept) -0.019 0.007 -0.033 -0.006 0.006
mr_plot(MRInputObject_prostate_lipid,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)